function gram_schmidt_qr()
% GRAM_SCHMIDT_QR Gram-Schmidt QR分解演示
% 
% 演示经典和改进的Gram-Schmidt正交化过程

fprintf('=== Gram-Schmidt QR分解 ===\n\n');

%% 1. 经典Gram-Schmidt过程
fprintf('1. 经典Gram-Schmidt过程\n');
fprintf('将矩阵A的列向量正交化，得到 A = Q*R\n\n');

% 创建测试矩阵
A = [1, 1, 0;
     1, 0, 1;
     0, 1, 1;
     0, 0, 1];

fprintf('原始矩阵 A =\n');
disp(A);
fprintf('矩阵维数: %d×%d\n', size(A, 1), size(A, 2));

% 经典Gram-Schmidt
[Q_cgs, R_cgs] = classical_gram_schmidt(A);

fprintf('经典Gram-Schmidt结果:\n');
fprintf('Q矩阵 =\n');
disp(Q_cgs);
fprintf('R矩阵 =\n');
disp(R_cgs);

% 验证正交性
QtQ = Q_cgs' * Q_cgs;
fprintf('Q^T*Q - I 的误差: %e\n', norm(QtQ - eye(size(QtQ)), 'fro'));

% 验证分解
QR_product = Q_cgs * R_cgs;
fprintf('Q*R - A 的误差: %e\n', norm(QR_product - A, 'fro'));

%% 2. 改进的Gram-Schmidt过程
fprintf('\n2. 改进的Gram-Schmidt过程\n');
fprintf('数值稳定性更好的Gram-Schmidt变形\n\n');

[Q_mgs, R_mgs] = modified_gram_schmidt(A);

fprintf('改进Gram-Schmidt结果:\n');
fprintf('Q矩阵 =\n');
disp(Q_mgs);
fprintf('R矩阵 =\n');
disp(R_mgs);

% 验证正交性
QtQ_mgs = Q_mgs' * Q_mgs;
fprintf('Q^T*Q - I 的误差: %e\n', norm(QtQ_mgs - eye(size(QtQ_mgs)), 'fro'));

% 验证分解
QR_mgs = Q_mgs * R_mgs;
fprintf('Q*R - A 的误差: %e\n', norm(QR_mgs - A, 'fro'));

%% 3. 数值稳定性比较
fprintf('\n3. 数值稳定性比较\n');
fprintf('使用病态矩阵测试两种方法的稳定性\n\n');

% 构造病态矩阵
epsilon = 1e-10;
A_ill = [1, 1, 1;
         epsilon, 0, 0;
         0, epsilon, 0;
         0, 0, epsilon];

fprintf('病态矩阵 A (ε = %e):\n', epsilon);
disp(A_ill);
fprintf('条件数: %.2e\n', cond(A_ill));

% 经典Gram-Schmidt
[Q_cgs_ill, R_cgs_ill] = classical_gram_schmidt(A_ill);
QtQ_cgs_ill = Q_cgs_ill' * Q_cgs_ill;
orthogonality_error_cgs = norm(QtQ_cgs_ill - eye(size(QtQ_cgs_ill)), 'fro');

% 改进Gram-Schmidt
[Q_mgs_ill, R_mgs_ill] = modified_gram_schmidt(A_ill);
QtQ_mgs_ill = Q_mgs_ill' * Q_mgs_ill;
orthogonality_error_mgs = norm(QtQ_mgs_ill - eye(size(QtQ_mgs_ill)), 'fro');

fprintf('正交性误差比较:\n');
fprintf('经典Gram-Schmidt: %e\n', orthogonality_error_cgs);
fprintf('改进Gram-Schmidt: %e\n', orthogonality_error_mgs);
fprintf('改进方法的优势: %.1f倍\n', orthogonality_error_cgs / orthogonality_error_mgs);

%% 4. 几何解释
fprintf('\n4. 几何解释\n');
fprintf('Gram-Schmidt过程的几何意义是逐步构造正交基\n\n');

% 使用简单的2D例子进行可视化
A_2d = [1, 1; 2, 0];
fprintf('2D示例矩阵:\n');
disp(A_2d);

[Q_2d, R_2d] = modified_gram_schmidt(A_2d);
fprintf('正交化后的Q矩阵:\n');
disp(Q_2d);

% 绘制几何解释
figure('Name', 'Gram-Schmidt几何解释');
subplot(1, 2, 1);
% 原始向量
quiver(0, 0, A_2d(1, 1), A_2d(2, 1), 0, 'r', 'LineWidth', 2);
hold on;
quiver(0, 0, A_2d(1, 2), A_2d(2, 2), 0, 'b', 'LineWidth', 2);
axis equal; grid on;
xlim([-0.5, 2.5]); ylim([-0.5, 2.5]);
title('原始向量');
legend('a₁', 'a₂', 'Location', 'best');

subplot(1, 2, 2);
% 正交化后的向量
quiver(0, 0, Q_2d(1, 1), Q_2d(2, 1), 0, 'r', 'LineWidth', 2);
hold on;
quiver(0, 0, Q_2d(1, 2), Q_2d(2, 2), 0, 'b', 'LineWidth', 2);
axis equal; grid on;
xlim([-0.5, 1.5]); ylim([-0.5, 1.5]);
title('正交化后的向量');
legend('q₁', 'q₂', 'Location', 'best');

%% 5. QR分解的唯一性
fprintf('\n5. QR分解的唯一性\n');
fprintf('当A列满秩且要求R的对角元为正时，QR分解唯一\n\n');

% 检查R的对角元
diag_R = diag(R_mgs);
fprintf('R的对角元: [%.6f, %.6f, %.6f]\n', diag_R);
fprintf('所有对角元为正: %s\n', mat2str(all(diag_R > 0)));

% 与MATLAB内置QR分解比较
[Q_matlab, R_matlab] = qr(A, 0);  % 经济型QR分解
fprintf('与MATLAB QR分解的差异:\n');
fprintf('||Q_mgs - Q_matlab||: %e\n', norm(Q_mgs - Q_matlab, 'fro'));
fprintf('||R_mgs - R_matlab||: %e\n', norm(R_mgs - R_matlab, 'fro'));

%% 6. 性能分析
fprintf('\n6. 性能分析\n');
fprintf('比较不同QR分解方法的计算时间\n\n');

sizes = [50, 100, 200, 400];
fprintf('矩阵规模    经典GS    改进GS    MATLAB QR\n');
fprintf('--------    -------   -------   ---------\n');

for i = 1:length(sizes)
    m = sizes(i);
    n = min(m, 50);  % 保持列数不太大
    A_test = randn(m, n);
    
    % 经典Gram-Schmidt
    tic;
    [Q_cgs_test, R_cgs_test] = classical_gram_schmidt(A_test);
    time_cgs = toc;
    
    % 改进Gram-Schmidt
    tic;
    [Q_mgs_test, R_mgs_test] = modified_gram_schmidt(A_test);
    time_mgs = toc;
    
    % MATLAB QR
    tic;
    [Q_matlab_test, R_matlab_test] = qr(A_test, 0);
    time_matlab = toc;
    
    fprintf('%6d      %5.3f s   %5.3f s   %7.3f s\n', ...
            m, time_cgs, time_mgs, time_matlab);
end

%% 7. 重正交化
fprintf('\n7. 重正交化\n');
fprintf('对于严重病态的矩阵，可能需要重正交化\n\n');

% 构造极度病态的矩阵
A_severe = [1, 1, 1;
            1e-15, 0, 0;
            0, 1e-15, 0;
            0, 0, 1e-15];

fprintf('极度病态矩阵条件数: %.2e\n', cond(A_severe));

% 标准改进Gram-Schmidt
[Q_std, R_std] = modified_gram_schmidt(A_severe);
orthogonality_std = norm(Q_std' * Q_std - eye(size(A_severe, 2)), 'fro');

% 重正交化Gram-Schmidt
[Q_reorth, R_reorth] = reorthogonalized_gram_schmidt(A_severe);
orthogonality_reorth = norm(Q_reorth' * Q_reorth - eye(size(A_severe, 2)), 'fro');

fprintf('正交性误差比较:\n');
fprintf('标准改进GS: %e\n', orthogonality_std);
fprintf('重正交化GS: %e\n', orthogonality_reorth);

%% 8. 应用：最小二乘求解
fprintf('\n8. 应用：最小二乘求解\n');
fprintf('使用QR分解求解最小二乘问题\n\n');

% 超定系统
m = 10; n = 4;
A_ls = randn(m, n);
x_true = randn(n, 1);
b_ls = A_ls * x_true + 0.1 * randn(m, 1);  % 添加噪声

fprintf('最小二乘问题: %d×%d 矩阵\n', m, n);

% 使用QR分解求解
[Q_ls, R_ls] = modified_gram_schmidt(A_ls);
x_qr = R_ls \ (Q_ls' * b_ls);

% 与正规方程比较
x_normal = (A_ls' * A_ls) \ (A_ls' * b_ls);

% 与MATLAB内置求解器比较
x_matlab_ls = A_ls \ b_ls;

fprintf('解的比较:\n');
fprintf('QR方法: [%.6f, %.6f, %.6f, %.6f]\n', x_qr);
fprintf('正规方程: [%.6f, %.6f, %.6f, %.6f]\n', x_normal);
fprintf('MATLAB: [%.6f, %.6f, %.6f, %.6f]\n', x_matlab_ls);

% 残差比较
residual_qr = norm(A_ls * x_qr - b_ls);
residual_normal = norm(A_ls * x_normal - b_ls);
residual_matlab = norm(A_ls * x_matlab_ls - b_ls);

fprintf('残差比较:\n');
fprintf('QR方法: %e\n', residual_qr);
fprintf('正规方程: %e\n', residual_normal);
fprintf('MATLAB: %e\n', residual_matlab);

end

function [Q, R] = classical_gram_schmidt(A)
% CLASSICAL_GRAM_SCHMIDT 经典Gram-Schmidt正交化
% 
% 输入: A - m×n 矩阵
% 输出: Q - m×n 正交矩阵, R - n×n 上三角矩阵

[m, n] = size(A);
Q = zeros(m, n);
R = zeros(n, n);

for j = 1:n
    v = A(:, j);
    
    % 减去之前所有正交向量的投影
    for i = 1:j-1
        R(i, j) = Q(:, i)' * A(:, j);
        v = v - R(i, j) * Q(:, i);
    end
    
    % 归一化
    R(j, j) = norm(v);
    if R(j, j) < eps
        error('矩阵A线性相关');
    end
    Q(:, j) = v / R(j, j);
end

end

function [Q, R] = modified_gram_schmidt(A)
% MODIFIED_GRAM_SCHMIDT 改进的Gram-Schmidt正交化
% 
% 输入: A - m×n 矩阵
% 输出: Q - m×n 正交矩阵, R - n×n 上三角矩阵

[m, n] = size(A);
Q = A;  % 复制A作为工作矩阵
R = zeros(n, n);

for j = 1:n
    % 归一化第j列
    R(j, j) = norm(Q(:, j));
    if R(j, j) < eps
        error('矩阵A线性相关');
    end
    Q(:, j) = Q(:, j) / R(j, j);
    
    % 从后续列中减去投影
    for k = j+1:n
        R(j, k) = Q(:, j)' * Q(:, k);
        Q(:, k) = Q(:, k) - R(j, k) * Q(:, j);
    end
end

end

function [Q, R] = reorthogonalized_gram_schmidt(A)
% REORTHOGONALIZED_GRAM_SCHMIDT 重正交化Gram-Schmidt
% 
% 对每个向量进行两次正交化以提高数值稳定性

[m, n] = size(A);
Q = A;
R = zeros(n, n);

for j = 1:n
    % 第一次正交化
    for i = 1:j-1
        r_ij = Q(:, i)' * Q(:, j);
        R(i, j) = R(i, j) + r_ij;
        Q(:, j) = Q(:, j) - r_ij * Q(:, i);
    end
    
    % 第二次正交化（重正交化）
    for i = 1:j-1
        r_ij = Q(:, i)' * Q(:, j);
        R(i, j) = R(i, j) + r_ij;
        Q(:, j) = Q(:, j) - r_ij * Q(:, i);
    end
    
    % 归一化
    R(j, j) = norm(Q(:, j));
    if R(j, j) < eps
        error('矩阵A线性相关');
    end
    Q(:, j) = Q(:, j) / R(j, j);
end

end